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ABSTRACT 


Three  approximation  methods  for  generating  outcomes  on 
Poisson  random  variables  are  discussed.   A  comparison  is  made 
to  determine  which  method  requires  the  least  computer  execu- 
tion time  and  to  determine  which  is  the  most  robust  approxima- 
tion.  Results  of  the  comparison  study  suggest  the  method  to 
choose  for  the  generating  procedure  depends  on  the  mean  value 
of  Poisson  random  variable  which  is  being  generated. 
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I.   INTRODUCTION 

It  is  frequently  desired  to  generate  Poisson  random 
variates  in  simulations.   There  are  standard  exact  methods 
for  doing  this;  the  problem  arises  when  a  computer  is  used 
to  generate  the  Poisson  random  number  which  has  a  large  mean. 
For  example,  generating  one  random  number  such  as  105  from  a 
Poisson  distribution  with  mean  100  needs  at  least  105  calls 
to  a  pseudo-random  number  (uniform  (0,1))  generator.   Compu- 
ter time  requirements  become  important  cost  factors  when 
considering  various  methods  for  generating  random  numbers. 

The  objective  of  this  paper  is  to  examine  several  ap- 
proximated ways  of  generating  Poisson  random  variables  and 
to  determine  the  method  which  gives  minimum  execution  time 
and  small  mean  absolute  deviation  according  to  the  Poisson 
mean  value.   The  mean  value  is  the  only  parameter  in  the 
distribution.   Comparison  statistics  to  determine  the  best 
approximation  to  the  Poisson  distribution  are  the  cumulative 
probability,  mean  absolute  deviation  and  Kolmogorov-Smirnov 
test.   Finally  a  composite  generating  procedure  according  to 
the  mean  value  is  suggested. 

The  following  notation  is  used  in  this  study. 

U.  denotes  a  uniform  (0,1)  random  variable; 

N   denotes  a  Poisson  random  variable  where  mean  is  m; 

Z   denotes  a  random  variable  from  a  unit  normal  distribution 


II.   GENERATION  OF  POISSON  DISTRIBUTED  VARIATES 

A.   THE  POISSON  DISTRIBUTION 

A  random  variable  N  with  integer  values  has  a  Poisson 

distribution  if 

-m  n 
prob  {N  =  m}  =  e  ?    n  =  0,1,2...   . 

In  order  to  generate  a  Poisson  random  number  N  from  a  Poisson 
distribution  with  mean  m,  the  following  algorithm  is  pre- 
sented.  It  is. the  standard  exact  method  for  generating  these 
deviates . 

Let  U.,  i  =  1,2,...,  be  independent  uniform  (0,1)  random 
variates .   The  Poisson  variate,  N,  is  computed  as: 


0   if  U.  <  e  m 
l  - 

N  = 

k  k+1 

k  if  n  U.  >  e   >  n  U. 

i=l  i=l 


where  N  is  non  distributed  as  a  Poisson  with  mean  m,  i.e., 
Prob  (N=n)  =  e~mmn/n!  . 

Equivalent ly,  since  the  logarithm  is  a  monotone  transforma- 
tion, we  have 

0   if  ln(U. )  <  ln(e~m)  =  -m 


N 


k         k 
k   if  ln(  31   U  )  =  I      ln(U  )  >  ln(e~m) 
i=l   1    1=1 

k+1  k+1 

=  -m  >   E   ln(U. )  =  In  (  n   U. ) . 
i-1  i=l 


Changing  the  signs  and  direction  of  the  inequalities  we 
have : 


0   if  -ln(U. )  >  m 
1   - 


N  = 


k+1 


k    if   £   -  ln(U.)  <  m  <   Z   -  ln(U. ) 
1=1       x      "1-1 

or  letting  E.  =  -  ln(U. ) 

°     l         l 

0    if  E.  >  m 
l  - 

k  k+1 

if  E  E.  <  m  <  E  E. 
•  t  i  -  .  ,  i 
i=l  i=l 

where  the  E.  are  exponentially  distributed  variates  with 
mean  1. 

If  n  multiplications  of  uniform  (0,1)  random  numbers  is 
strictly  greater  than  e~   and  if  n+1  multiplication  of  uniform 
(0,1)  random  number  is  equal  and  less  than  e    then  n  is  the 
Poisson  random  number.   Generally  generating  one  random  num- 
ber from  Poisson  process  with  parameter  m  requires  on  the 
average  m+1  uniform  (0,1)  random  numbers.   This  is  because 
the  number  generated  is  n+1.   When  m  is  large  it  is  clear 
that  generating  Poisson  random  numbers  with  the  above  method, 
although  it  is  exact  takes  a  lot  of  computer  time  and  this 
method  may  be  uneconomical.   In  addition,  the  large  number 
of  multiplications  can  produce  serious  precision  problems 
on  a  digital  computer. 
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B.   APPROXIMATIONS  FOR  POISSON  VARIATES 
1.   Normal  Approximation 

In  a  Poisson  process  with  parameter  X  it  is  necessary 
to  generate  random  variables  from  the  Poisson  distribution 
with  parameter  (mean)  m.   Now  look  at  counts  in  (03x),  where 
x  satisfies  Xx  =  m.   The  central  limit  theorem  says  that  as 
m  goes  to  positive  infinity,  (or  when  x  gives  the  infinity  in 
Poisson  process  with  fixed  X),  then  N,  which  has  mean  m  and 
variance  m,  is  such  that  N+^-m/m  is  approximately  distributed 
as  a  unit  normal  random  variable.   Denote  a  random  variable 
from  a  unit  normal  distribution  by  Z .   So  N  is  distributed 
approximately  as  m^Z  +  m.   In  order  to  generate  Poisson  ran- 
dom numbers  from  the  normal  distribution,  first  generate  Z; 
then  let 

0  if  m^Z  +  m  -  0.5  <  1 


N  = 

| m^Z  +  m  -  0 . 5 1        otherwise 

where  | a|  denotes  the  greatest  integer  less  than  or  equal  to 
a;  also  known  as  the  "floor"  of  a.   N  is  then  the  approximated 
Poisson  random  variate. 

2 .   Square  Root  -Transformation  of  Poisson  Distribution 
If  N  is  a  Poisson  random  variable  with  mean  m  then 
Y  =  /N  +  3/8  is  approximately  distributed  as  a  normal  dis- 
tribution  with  mean  m   and  variance  \.      This  result  is  due  to 
Bartlett  [Ref.  9], 
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This  method  is  derived  as  follows:   let  Y  =  /N  +  C   ~  N(y,a2) 
where  C  is  a  non-negative  constant.   Let  t  =  N  -  m  and 
m1  =  m  +  C.   Define  coefficients  for  s  =  1,2,3,... by 

=  (-D: 


a    ^_^s+l  l'(-l)-(-3)---(-2s  +  3) 


2s  •  s ! 


Then  for  any  t  >  -  m'  we  have  a  Taylor  series  expansion, 


Y  =  /m*"  {1+A1  ^r  -  A0(^r)2  +.--  +  (-l)sA  n(^-r)s+1} 
1  m'     2  m!        v  '      s-1  m' ' 

If  t  >  0,  we  see  at  once  |R  |  <  A  t  /(m')   5  converges  and 

5        S 

is  bounded  (F.  J.  Anscombe  (4)).   We  note  now  that  the  moments 
of  t  are  ]i1    =  0,  y2  =  m>  y  3  =  m,  yi,  =  3m2  +  m,...,  which  give 


so  that  when  C  =  3/8  Var  (Y)  -  %    (1  +  l/l6m2 ) .   Also 


E(Y)  ~  ^+C \-  +    24C  "^  . 


8m3*    128m3/2 

Let  XNR  =  /N  +  3/8.   Then  XNR  is  approximately  normally  dis- 
tributed with  mean  /m  and  variance  \. 


Z 

XNR 

= 

XNR 

-    vm 

+    3/8 

/m   + 

3/8    , 

/N" 
thus    se 

+ 
it 

3/8 

/m   + 

3/8    , 
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0  if  (|  +  /m+3/b)2  -  3/8  <  1 

N  -  " 

(|  +  /m+3/8)2-  3/8 I  otherwise 

N  is  then  the  approximated  Poisson  variate. 

We  now  need  to  calculate  the  probability  distribution 
of  N  obtained  in  this  way  from  the  square  root  transformation, 
We  want  the  probability  that 

n  -  1  +  3/8  <  XNR2  <  n  +  3/8 
if  we  divide  by  the  variance  \   we  get 

4(n  -  1  +  3/8)  <  XNR2  4  <  4(n  +  3/8). 
Note  that  4*XNR2  is  distributed  with  a  non-central  x2  distri- 
bution.  The  non-central  Xi  density  is 

e-to+y2/c2)        -  y ^  H/r 

f    (x)    =   [e      a        +   e    a         ]. 

x  2/x      /2i, 

Thus    if  y=m^3    o2=(h)2 s    then 

f    (x)    =    e^<^«    [e-(^A)^  +   e(^A), 
X  2/x      /2li 


and 


N+3/8  N+3/8-1 

J(:i)2 
prob(x=N)  =  \"      f  (x)dx  -   ;        f  (x)dx 


yU)    fx(X)dx  -   0( 


This  allows  us  to  evaluate  directly  how  well  the 
distribution  of  N  approximates  the  distribution  of  a  Poisson 
variate  with  parameter  m. 

13 


Note  that  since  in  the  LLRANDOM  package  it  takes  the 
same  amount  of  time  to  generate  5  uniform  random  variables 
as  it  takes  to  generate  a  normal  random  variable,  the  pro- 
cedures will  be  competitive  timewise  once  m  is  much  greater 
than  5 • 

3.   Cube  Root  Transformation  of  Pois-son  Distribution 

If  N  is  a  Poisson  random  variable  with  mean  m  then 


Y  =  i^N- 1/24  is  approximately  distributed  as  a  normal  distri- 
but ion  with  mean  m    and  variance  1/9  3Vm  when  N  >  1.   This 
comes  from  the  following  derivation  which  is  essentially  the 

same  procedure  as  for  the  square  root  transformation.   Sup- 

3 

pose  Y  =  /N+C  is  distributed  as  N(y,a2).   Let  t  =  N  -  m  and 

m'  =  m  +  C.   Define  coefficients  for  s  =  1,233»..-  by 

_   _  (-l)s+1  i.(-2)(-5)(-8)-»y4-3s) 

a   — < . 

s  0s   , 

3   s! 

For  any  t  >  -  m'  we  have  the  Taylor  series  expansion 

if  t  >  0,  we  see  at  once  that  I R  \£.   t  /(m')s~   ^  converges. 

'  s  "^s  s 

Therefore , 

f(s)  (i  +e  5T  )  t  B 
Rs  "  si       (S>  °  S  e  <  1 

|t|  <  m* 
Rs(m')-1/3  =  (1  +  £   )1/3  -  U  +  ax  ^  -..  +  (-l)s  a^)8"1  ) 


Rs(m')"1/3    =0 

=   £  (-1)    a. (— r)   converges  and  is 

tS  1-8  X  m 


bounded. 


We  note  now  that  the  moments  of  t  are  y1=0,  \i2=mi 
u3=m,  u  =3m  +m,...,  giving  us 


1  1 


and 


E(Y)  =  3v^+C  -  J^  -±j~   +.  .. 

in 

Var(Y)  =    ^/m"5")2  {  j'j^jz    ™  - 


1  2m2+m 
81  (m1)* 


^v^m 


1  (      y   .,   1   >, 


3  .>ct   27 


/m1 


16T 


If  c  = 


=  24" 


then 


Var(Y)  = 


93/m 


Let 


YNR  =  3/N  -  1/24    N  >  1 


„  _  YNR  -  Vm  +  1/24 

36^ 


YNR  =  3/m  +  1/24  +(l/36 /m)  Z 


3 /N- 1/2 4  =  3/m  +  1/24  +  (l/36/m)  Z 


Thus  set 


N  = 


|  ( 3  /m+1/24"  +  3  6  ^  Z  )  3  +  ^ 


if  (3/m+l/24+i  6/m"  Z)3+^  <  1 


otherwise 
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N  is  now  the  approximated  Poisson  variate. 
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III.   EVALUATION  OF  THE  METHODS 

Generally,  generating  Polsson  random  number  from  the 
exact  method  is  known  to  take  a  long  execution  time,  since 
one  generated  random  number  "N"  requires  on  the  average  N 
multiplications  of  uniform  (0,1)  random  numbers.   Therefore, 
generating  Poisson  variables  with  the  exact  method  (which 
gives  the  best  accuracy)  is  good  for  a  small  mean,  m,  while 
the  approximation  methods,  which  take  shorter  execution  time 
but  with  less  accuracy,  should  be  used  for  large  m.   Here  we 
need  a  trade-off  between  execution  time  and  accuracy  to 
choose  the  generation  method  according  to  the  mean  value  of 
the  Poisson  distribution. 

The  comparison  statistics  show  how  closely  the  methods 
approximate  the  original  Poisson  distribution.   In  the 
Kolmogorov-Smirnov  test,  all  approximations  are  accepted  at 
significance  level  a  =  .05.   Prom  the  comparison  of  the  em- 
pirical probability  distributions  and  the  mean  absolute 
deviations  from  the  exact  distribution,  the  optimal  genera- 
tion procedure  based  on  the  mean  value,  m,  is  as  follows: 

Method  Mean  (mj 

1.  exact  Poisson  distribution  if  0  <  m  <  20 

2.  square  root  transformation  -if  20  <  m  <  100 

3.  normal  approximation  if  m  >  100. 

The  cube  root  transformation  should  not  be  adopted  be- 
cause it  is  far  less  accurate  compared  to  the  square  root 
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transformation  and  the  normal  approximation.   The  following 
tables  and  figures  show  the  comparison  statistics  of  the 
three  approximations  versus  the  exact  Poisson  distribution. 

Table  I  assesses  the  accuracy  of  the  normal,  square  root, 
and  cube  root  approximations  to  the  exact  Poisson  distribu- 
tion at  selected  points.   Figure  1  illustrates  graphs  of  the 
empirical  cumulative  distribution  functions  of  the  three 
approximations  to  that  of  an  exact  Poisson  distribution  with 
mean  10.   Figure  2  shows  graphs  of  the  mean  absolute  devia- 
tions of  the  three  approximations  from  the  exact  Poisson 
cumulative  distribution  function  for  various  mean  values. 
The  mean  absolute  deviation  is  defined  as: 

<ra  jx  <pi  -v2> 

where 

P!  is  the  CDF  of  the  approximating  distribution; 

P.  is  the  CDF  of  the  exact  Poisson  distribution:  and 
l 

k   is  the  sample  size. 

Table  II  is  a  summary  of  a  comparison  of  the  empirical 
distributions  produced  by  the  three  approximation  methods  to 
the  exact  Poisson  distribution  by  means  of  the  one  sample 
Kolmogorov-Smirnov  test.   The  null  hypothesis  is  that  the 
approximated  distributions  are  Poisson  against  the  alterna- 
tive that  they  are  not  Poisson. 

Lastly,  Table  III  and  Figure  3  analyze  the  sensitivity  of 
the  square  root  transformation  to  the  constant  C.   In  the 
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derivation  of  the  square  root  transformation,  C  was  chosen 
as  3/8.   Different  constants  were  used  in  order  to  find  the 
most  robust  constant  to  use,  i.e.,  the  constant  which  yielded 
the  smallest  mean  absolute  deviation  from  the  exact  Poisson. 
The  value  C  =  13/18  was  found  to  be  the  most  robust.   Note 
that  this  value  of  C  was  used  in  the  approximation  when 
making  the  comparisons  with  the  other  methods. 
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Table  I.   Accuracy  of  the  Normal,  Square  Root,  and  Cube  Root 
Approximations  to  the  Exact  Poisson  Distribution. 


Poisson 
mean,  m 

P<Z(m)  <  n> 

Observed 
Value,  n 

Exact 

Normal 
Approx 

Square  Root 
Transformation 
C=13/18 

Cube  Root 

Transformation 

I   5 

2' 

3 
6 

0.08422 

0.14037 
0.14622 

0.07301 
0.11939 
0.16036 

0.07719 
0.14051 
0.14614 

0.10491 
0.18491 
0.16681 

15 

10 
14 
18 

0.04861 
0.10244 
0.07062 

0.04485 
0.09937 
0.07622 

0.04704 
0.10288 
0.07073 

0.09045 
0.05512 
0.06348 

20 

14 

19 
22 

0.03874 
0.08884 
0.07692 

0.03633 
0.08683 
0.08050 

0.03746 
0.08891 
0.07672 

0.08819 
0.08726 
0.05755 

40 

36 

39 
42 

0.05394 
0.06296 
O.05850 

0.05161 
0.06223 
0.05995 

0.05409 
0.06307 
0.05841 

0.02552 
0.06336 
0.05675 

61 

55 
59 

0.03960 
0.05019 

0.03910 
0.04999 

0.03963 
0.05020 

0.00662 
0.04973 

;  82 

82 
87 
89 

0.04407 
0.03686 
0.03165 

0.04402 
0.03649 
0.03243 

0.04402 
0.03683 
0.03149 

0.03289 
0.03801 
0.03991     I 

i  90 

94 
97 

0.03775 
0.03H2 

0.03805 
0.03182 

0.03759 
0.03098 

0.0255 
0.0435 

100 

93 
99 

0.03223 
0.03997 

0.03241 
0.03980 

0.03217 
0.03994 

o.onio 

0.03605 

120 

112 
116 

0.02881 
0.03477 

0.02847 
0.03438 

0.02866 
!   0.03462 

0.03627 
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1   c 

Mean  Absolute  Devi 

/l   k 

ationl/—  E  (PJ  ■ 
Kk-l1=1   i 

-  v2 

Polsson  Mean  =  20 

Polsson  Mean  = 

80 

3/8 

0.019 

0.009 

|  5/9 

0.016 

0.006 

j  7/12 

0.014 

0.005 

11/18 

0.011 

0.005 

23/36 

0.009 

0.004 

8/12 

0.007 

.  0.003 

•13/18 

0.006 

0.002 

|  9/12 

0.007 

0.003 

5/6 

0.011 

0.004 

*  optimal  constant  for  C 


Table  III. 


Sensitivity  of 'the  Square  Root  Transformation  to  the  Constant 
C  for  Poisson  Distributions  with  Means  m  =  20  and  m  =  80 . 
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Figure  3.  Plot  of  Sensitivity  Test  Given  in  Table  III. 
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IV.   THE  METHODS  OF  DIETER  AND  AHRENS 

During  the  writing  of  this  thesis,  a  pre-publication 
chapter  from  a  book  by  U.  Dieter  and  J.  H.  Ahrens  was  re- 
ceived.  This  chapter  deals  with  the  generation  of  non-uni- 
form variates  and  contains  a  section  on  the  generation  of 
Poissan  variates.   Some  of  the  methods  which'  they  described 
will  be  outlined  here.   They  have  not  been  programmed  or 
tested  against  the  composite  method  of  this  thesis. 

A.   UNIT  MEAN  METHODS 

It  is  well-known  that  the  sum  of  independent  Poisson 
deviates  with  mean  y1,y2,,j  is  Poisson  distributed  with  mean 
Ui+u2+...   Hence  the  following  algorithm  may  be  considered. 

1.  Split  the  mean  m  of  the  Poisson  distribution  into 
m'-<-[m]  and  f«-m-m'<l.   (Hence  [m]  means  take  the  integer  part 
of  m  and  f  is  then  the  fractional  part.) 

2.  Take  m'  samples  from  the  Poisson  distribution  of 
mean  1  using  some  suitable  method.   Let  s  be  the  sum  of  these 
samples . 

3.  Take  a  further  sample  n  from  the  Poisson  (1)  distri- 
bution.  If  n  =  0  deliver  k-<-s .   Otherwise,  generate  U1}U2...U 
and  count  the  number,  j,  of  U.  for  which  U.  <  f,  and  deliver 
k^s+j . 

The  validity  of  step  3  is  proved  as  follows:  the  probability 
of  observing  j  samples  U.  less  than  or  equal  to  f  (once  n  -  j 
samples  greater  than  f)  is 
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n 


(  *  )  f J  (l-f)n~J. 

Since  n  follows  a  Poisson  distribution  with  mean  1  this 
yields 

°°   —  1  °°      —1 

P,  =  z     ^—r   (  ?  )  f J  (l-f)n-J  =  z  ,t/e  .,,  fJ'(i-f)n_J 
J    n=j   nI    J  n=j  TTU^JT 

=  g^fj   "    (i-f)n"J  =  e"1^'   1-f  =  e"ffJ' 

j'1    n-j=o    (n  d  l        y-   e  J [ 

The  unit-mean  method  becomes  a  suitable  algorithm  for 
the  case  when  the  exact  method  should  not  be  used  for  gen- 
erating Poisson  random  variatesa  that  is,  when  the  mean  is 
large  (m  >  200).   This  method  will  then  become  slightly  faster 
than  the  exact  method,  but  still  needs  to  call  too  many  uni- 
form deviates  since  it  is  basically  the  same  as  the  exact  method 
except  for  splitting  the  mean,  compared  to  the  square  root 
transformation  approximation  method,  the  unit  mean  method  may 
be  uneconomical. 

B.   THE  CENTER  TRIANGLE  METHOD 

This  method  is  suitable  in  the  case  of  variable  means  m 
of  the  Poisson  distribution.   The  mean  is  split  into  m*  =  [m] , 
the  integer  part  of  m,  and  the  remaining  fraction  0<f=m-m'<l. 
The  algorithm  uses  a  mixture  of  two  distributions.   The  first 
distribution  has  an  easy  sampling  method  and  is  employed  most 
of  the  time.   It  is  based  on  the  following  isosceles  triangle. 
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J  b    _b 
\  s    s2 

V  0    el: 


g(x)  =   /  —  -  —  |x-m' |  if  m'-s  <  x  <  m'+s 
(B.l) 

.sewhere 


The  constant  s  is  chosen  as 
(B.2)      s  =  2.2160358671665  /m7"  -  0.78125. 

The  area  under  g(x)  is  b.   If  b  is  small  enough  then  the 
inequalities 


.-J1 


+1 


(B.3)     t,  =  t\         g(x)dx  <  P.  (i  =  0,1,2,...) 


may  hold  for  all  i.   The  maximum  b  for  which  this  is  true  was 
found  numerically  for  all  4  <  m  <  200.   The  choice  of  the 
constant  a  =  2.2160358671665  in  (B.2)  is  well  motivated.   The 
vertices  of  the  maximum  isosceles  triangle  that  can  be  placed 
wholly  inside  the  standard  normal  curve  are  situated  at  (±a,0) 
and  (0,1//2tF).   As  m*  goes  to  Infinity,  b    approaches  the 

ffldX 

area  of  this  triangle  which  is  .88^407040229876.   The  correc- 
tion -0.78125  =  -25/32  in  (B.2)  improves  the  value  of  b  ^v 

max 

for  small  m.   [7,  page  11-15,  table.]   Here  the  variable  mean 
m  is  split  into  m  =  m'tf  as  in  the  unit  mean  method.   If 
m  <  8  the  exact  method  is  applied.   If  m  >  9j  the  triangular 
method  is  used  with  the  triangular  probabilities  t.  =  27/32. 
This  center  triangle  method  needs  a  long  program  including 
several  tables  which  requires  a  large  memory  space  and  suf- 
fers from  the  roundoff  errors '  of  b . 
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This  method  proves  a  very  efficient  algorithm  based  on 
a  simple  idea.   However,  it  is  somewhat  messy  to  describe 
the  program  and  has  little  mathematical  appeal  compared  to 
the  square  root  transformation  method.   This  method  also 
needs  a  long  execution  time  as  does  the  exact  method.   For 
example,  to  compute  one  Poisson  variate  from  a  distribution 
with  mean  m  =  100,  the  program  would  require  1886  ys  [7, 
page  11-19] . 

C.   THE  GAMMA  METHOD 

In  order  to  obtain  a  sample  k  from  the  Poisson  distribu- 
tion with  variable  mean  m,  select  a  positive  integer  or 
(typically  n  is  a  little  smaller  than  m) .   Then,  take  a  sam- 
ple x  from  a  Gamma  distribution  with  parameter  n. 

Case  (1)   if  x  >  m,  return  a  sample  k  from  the  binomial 
distribution  with  parameters  n-1,  m/x. 

Case  (2)  if  x  <  m,  take  a  sample  j  from  the  Poisson  dis- 
tribution of  mean  m-x  and  return  k<-n+j  . 

The  sample  x  simulates  the  n-th  event  (arrival)  in  a 
Poisson  process  of  rate  1.   If  x  >  m,  then  there  are  n-1 
arrivals  in  the  interval  (0,x),  and  each  of  these  has  a 
probability  of  m/x  of  being  below  m  (Case  (1)).   If  x  <  m, 
then  the  n  simulated  arrivals  are  all  before  m  and  the  sample 
j  indicates  the  additional  events  between  x  and  m  (Case  (2)). 

A  formal  proof  of  the  procedure  runs  as  follows .   In  the 
first  case  one  has  m  <  x  <  °°  and 
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~x   n—X 
Pk  =    (   k  }  (x}   (1"x}    ^     r(n)    dx 


is  the  probability  of  obtaining  k  from  the  binomial  distribu- 
tion (n-1,  m/x)  summed  over  the  Gamma  (n)-distributed  values 
of  x  above  m.   The  expression  transforms  into 


/(n-1)!     mk(x-m)n"1" 
k'(n-l-k)!       n-1 


i-1-k    -x   n-1 

Pk  =  J   k'(n-l-k)!   "  'V-i         (n-1)  I   ( ' K 
m 


s 


m      I   /    Nn-l-k   -x, 
kl(n-l-k)!  -"   (x~m)       e   dx 


m 


k   -m         ,      , 
me     |   ,n-l-K   -t 

t      e   dt 


kl(n-l-k) I ^ 

"  o 

where  t  =  x  -  m.   Notice  that  T(n-k)  is  (n-l-k) ! .   Hence 

k   -m 
p   =  m   e 
k     k! 

as  required.   In  the  second  case  one  has  0  <  x  <  m,  and 

f%  /    sk-n   -(m-x)    -x   n-1 
P   -  I   (m-x)     ej      e    x     , 

*k  ~J        (l^nTi        r(n)    ax 

is  the  probability  of  obtaining  j  =  k-n  from  the  Poisson 
distribution  of  mean  m-x  summed  over. the  Gamma  (^-distri- 
buted values  of  x  between  0  and  m.   This  expression  trans- 
forms into 
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3  e        f    n-1  /    vk-n  , 

i  =  t\ vi — 7 T\-r     I   x     (m-x)     dx. 

k    (k-n) !  (n-1) !  «l 


k    V,K-n;  ;  (.n-ij  :  , 

~o 

Introducing  t  =  x/m  yields 

1 


-m  n-1  k-n    „,     ,      , 
P  =  e   m    m    m   |   n-1  (1_t)k-n 

k     (n-1) ! (k-n) !    " 


S 


o 
The  integral  is  a  Beta  function  3(n,k-n+l) 

p   -    e"m  mk      r(n)  r(k-n+l)    e"m  mk 
'-   k  "  (n-1)!  (k-n)!       r(k+l)        k! 

as  before.    The  following  algorithm  is  considered. 

1.  Initialize  k-<-0 ,   w*-m 

2.  If  w  >  C  (C  =  24  was  used  as  the  mean  cut-off  point) 
go  to  6 

3.  (Start  Case  (2)).   Set  p*-l  and  calculate  b+-e~w 

4.  Generate  U  and  set  p«-pU.   If  p<b  deliver  k 

5.  Increase  k-*-k+l  and  go  to  4 

6.  Set  n<-[dwj  where  d  =  7/8.   Take  a  sample  x  from  the 
Gamma-(n)  distribution.   If  x>w  go  to  8. 

7.  Set  k  k+w,  w  w-x  and  go  to  2 

8.  (Start  Case  (1)).   Set  p^w/x 

9.  Generate  U.   If  U<p  increase  k-«-k+l 

10.  Set  n^-n-1.   If  n>l  go  to  9 

11.  Deliver  k. 

The  performance  of  the  algorithm  depends  on  the  cut-off 
point  C  in  Step  2.   Step  3»    Step  4  and  Step  5  are  exactly  the 
same  as  the  steps  of  exact  methods.   This  gamma  method  has  a 
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complex  sub  algorithm  for  the  binomial  distribution  in  Steps 
8-10  and  requires  memory  space  for  the  gamma  distribution. 
This  method  may  be  efficient  in  cases  with  extremely  large 
means,  .say  m  =  1000,  but  this  method  requires  a  rather  long 
execution  time  as  does  the  exact  method.   For  example,  to 
generate  a  sample  from  a  Poisson  distribution  with  mean  100 
takes  1175  us,  see  [7,  page  11-23]. 
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